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Multistable Bursting 

Neurons that support bursting dynamics are a common feature of neural sys- 
tems. Due to their prevalence, great effort has been devoted to understand- 
ing the mechanisms underlying bursting and information processing capabilities 
that bursting dynamics afford. In this paper, we provide a link between neu- 
ronal bursting and information storage. Namely, we show that the mechanism 
implicit to bursting in certain neurons may allow near instantaneous modifi- 
cations of activity state that lasts indefinitely following sensory perturbation. 
Thus the intrinsic, extra-synaptic state of these neurons can serve as a memory 
of a sensory event. 



I. INTRODUCTION 

Bursting is a dynamic state characterized by alternating periods of fast oscillatory behav- 
ior and quasi-steady-state activity. Nerve cells commonly exhibit autonomous or induced 
bursting by firing discrete groups of action potentials in time. Autonomously bursting neu- 
rons are found in a variety of neural systems, from the mammalian cortex- and brainstem^ 
to identified invertebrate neurons^. 

Multirhythmicity in a dynamical system is a specific type of multistability which describes 
the coexistence of two or more oscillatory attractors under a fixed parameter set. Multirhyth- 
micity has been shown to occur in vertebrate motor neurons^, invertebrate interneurons^, 
and in small networks of coupled invertebrate neurons^. Additionally, multirhythmicity 
has been demonstrated in models of intracellular calcium oscillations^ and coupled genetic 
oscillators^. 

Multistable systems can act as switches in response to an external input. The feasi- 
bility of multistability as an information storage and processing mechanism in neural sys- 
tems has been widely discussed in terms of neural recurrent loops and delayed feedback 
mechanisms^ - —. Theoretical studies have shown multirhythmic bursting behavior in a 
number of single neuron models^ 1 ^ as well as in a model two-cell inhibitory (half-center 
oscillator) network—. In biological neurons, it is possible these dynamics are employed as a 
short-term memory. In this report provide a general explanation for the existence of multi- 
rhythmic bursting in previous studies^ 1 ^ 1 ^ and the characteristics of a bursting neuron that 
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allow multirhythmic dynamics. Additionally, we provide experimental evidence suggesting 
the existence of this behavior in several identified bursting neurons of the aquatic mullusc 
Aplysia Californica. 

II. A SIMPLE PARABOLIC BURSTING MODEL 

Dynamical bursting systems are a subset of the singularly perturbed (SP) class of differ- 
ential equations, 

x = f(x,y), (1) 
y = eg(x,y),xeR m ,yeR n , (2) 

where < e is a small parameter. Using singular perturbation methods^ 1 ^, the dynamics 
of bursting models can be explored by decomposing the full system into fast and slow- 
subsystems: Eq. [1] and Eq. [2j respectively. The slow-subsystem can act independently^, be 
affected synaptically 17 , or interact locally with the spiking fast-subsystem*^ - — to produce 
alternating periods of spiking and silence in time. To examine the dynamical mechanism 
implicit to a bursting behavior, y is treated as a bifurcation parameter of the fast-subsystem. 
This is formally correct when e = and Eq. [2] degenerates into an algebraic equation, but 
the assumption is reasonable when there is large time separation between fast and slow 
dynamics. 

Using this technique, all autonomously bursting single neuron models displaying mul- 
tirhythmic bursting in the literature^ 1 ^ are topologically classified as the circle/circle 
type22i22; that is, their fast-subsystem is driven back and forth across a saddle node on 
invariant circle (SNIC) bifurcation to produce alternating spiking and silent states. These 
models can be reduced to a form that supports a topological normal SNIC to and from the 
active phased This system is, 

v = I + v 2 + Ul , (3) 
iii = —au 2 , (4) 
u 2 = -fi{u 2 - ut), (5) 

where J is a constant current, a and (3 are small positive constants. Trajectories are reset 
after a voltage spike by v — v c , v <— v T and (ui,u 2 ) {u\ + di,u 2 + d 2 ) where v c — v r , d\ 
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and g?2 are discrete shifts in variables that account for the effects of a spike. Thus, the slow- 
subsystem defined by Eq. H] and Eq. [51 is a damped linear oscillator defining a stable focus 
when (3 < 4a or node when (3 > 4a. Under the parameter sets used here, the slow-subsystem 
is a focus. Eq. [3] is the fast-subsystem and is a quadratic integrate and fire neuron. 



Eqs. (O HI |5]) is a SP system for small a and 0. When u% is used as a bifurcation parameter 
of the singular system, Eq. [3J, a saddle node bifurcation occurs when u\ = u sn = —I. When 
U\ > u sn , v — » oo like tan(t) (see Appendix IA II) . In the full system, when u\ < u SQ 
trajectories slowly converge on the equilibrium point of the slow focus. When u\ > u sn , 
trajectories of the slow subsystem are interrupted by spiking events when v goes to v c 
and (^1,^2) are discretely modified. The existence of a bursting solution is reliant on the 
interaction between spiking events and a flow of the slow focus - neither activity type can 
exist indefinitely when isolated. A necessary condition for the existence of a limit cycle 
that represents bursting or tonic spiking is that the equilibrium point of the slow-subsystem 
(ttijUjj) m ust satisfy u\ > u sn . 



A general aspect of SP systems of the form Eqs. (jH [2]) is that if M = {(x, y)\f(x, y) = 
0} is a stable equilibrium manifold of the fast subsystem on which D x f is non-singular, 
trajectories on M follow the reduced field, y = g(h(y),y), where h(y) is a function satisfying 
f{h(y),y) = 0^. With this in mind, consider a m + n dimensional circle/circle bursting 
system (e.g. Eqs. (El SI E})- Let <f> t be a flow of the full system. On M, h(y) is dependent 
on y and static parameters. With e = 0, <p t (x,y) = (j) t (h(y), y) since the times-scales of the 
singular and slow-subsystem are assumed to have infinite separation. Hence, if local cross 
section II C R m+n of dimension m + n — 1 is everywhere transverse to <f) t , the condition 
(fi T (h(yo), yo) = (h(yo),yo) on IT shows a r-periodic closed orbit. This condition reduces to 
<p T (yo) = yo and fixed points of a n — 1 dimensional map generated from an-1 dimensional 
section on M show closed orbits in R, m+n (figure [Q. In the case of Eqs. (JBl HI E])), n = 2. 
Therefore a one- dimensional return map equivalent to the full system can be created using 
a one-dimensional section S. 
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FIG. 1. A generic circle/circle bursting system projected into R 1+2 . The dynamics of x are much 
faster than those of y. At some points, trajectories are smashed onto the fast equilibrium manifold 
allowing return maps created from Poincare sections of dimension n — 1 to provide a complete 
dynamical description of the continuous system (see text for details). 

III. BURSTING SOLUTIONS RESULT FROM COUNTERACTING 
DYNAMICS 

Figure |2] indicates that discrete spiking events on a circle/circle bursting solution have 
a contraction-balancing action in phase space; they periodically force the stable focus, re- 
turning it to its initial condition after a single period, transforming the focal trajectory into 
a limit cycle. When trajectories containing different numbers of spiking events return to 
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their initial condition after a full revolution, multiple coexisting limit cycles are formed and 
multirhythmic bursting is achieved. 

Consider multirhythmic bursting produced by Eqs. (El HI E]) in figure EJ Let M c represent 
the equilibrium manifold of Eq. [3] parameterized by U\. We recall from section HT1 that a 
n — 1 = 1 dimensional section is needed to form a return map for the system defined by 
Eqs. ([31II1 [5]) so long as that section is on M c . Therefore we can define this section as the line 
of saddle-nodes that divides the silent and spiking regimes of phase space, S c = {(«i,m 2 ) G 
R | Mi = u sn }. Consequently, a one dimensional recurrance map, P : u 2 — > u 2 is formed by 
intersections with the Poincare section^, 

XL = {(ui,u 2 ) G R | Mi G X c ,«'i < 0}. (6) 

To further explain how bursting and multirhythmicity arises in circle/circle systems, we 
divide the full return map P into two components, G and H in terms of two Poincare 
sections, £_, and an additional section, 

S + = {(ui,u 2 ) G R | u\ G E c , Ui > 0}. (7) 

G is a mapping from some intial condition, 1X2(0) G X_ to a point on £+. H is then a 
mapping from G(w 2 (0)) back to E_ (figure[2]). Thus, over the course of a periodic trajectory, 
G accounts fo dynamics in the silent (contractive) region to the left of E c and H accounts for 
dynamics in the spiking (expansive) region to the right of E c . The functional composition 
of these maps, 

P(u 2 )=H(G(u 2 )), (8) 

is the discrete time Poincare recurrence map as previously defined. Using G and H, distance 
metrics, 

C I (u 2 ) = u 2 -G(u 2 ), (9) 

C s (u 2 ) = G(u 2 )-H(u 2 ), (10) 

C T (u 2 ) = C s (u 2 )-C r (u 2 ), (11) 

can be compared to describe contraction of trajectories in the resting, spiking, and combined 
regions of state space, respectively (see figure [2}. For points u 2 G £_ that C T {u 2 ) = 0, the 
contraction of the resting portion of the trajectory is balanced by the net expansion of the 
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FIG. 2. Multirhythmic bursting produced by Eqs. ([31 [H [5|). (a) Voltage traces showing a single 
period of three coexisting, limit cycles 10, 11, and 12 spikes. Note that the period of coexisting 
bursting solutions is different as indicated by three endpoint markers, (b) The corresponding limit 
cycles projected into the plane of the slow variables (thick) with the averaged solution shown in 
grey. Dotted grey lines are the nullclines of the averaged bursting system (see Eq. [13] and Eq. fH|) . 
Contraction metrics, C r , C s , and C T (Eqs. (f9l PTOl [TT]) are demonstrated on a portion of a trajectory 
of the full system spiraling outward (thin). The dashed line (S + ) and dot-dashed line (£_) are the 
Poincare sections used to define G and H (Eqs. ([H [7]). (c) C r , C s , and C T for the multirhythmic 
system in (a). When the contraction of the slow focus is balanced by expansive spiking events, 
C T = and there is a fixed point of the map. Contraction metrics for the averaged system are thin 
solid lines. Parameters used here are I = 0.5 (I = 1.2 for averaged system), a = 0.2, f3 = 0.05, 
di = 0.4,d 2 = 0.6, v c = 10, v v = -1. 



spiking portion. Hence, u\ are fixed points of P and therefore show closed orbits in the full 
system. 

Contraction and expansion of Eqs. ([3j HI [5]) are clarified by examining the averaged slow- 
subsystem. This uses a near-identity change of variables to account for the time-averaged 
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,ir z ) 

/ g(v(t,z),z)dt, (12) 
Jo 



effect of fast spiking when u\ > u SQ . Let ip(t,u) be the limit cycle of the fast-subsystem 
defining T-periodic spiking, u = eg(<p(t, u),u) is the periodically forced slow-subsystem. By 
the Pontryagin-Rodygin theory^, the averaged slow-subsystem for u\ > u sn is defined as, 

rT[z) 

T(z) ' 

where z = u + 0(e). So, for our simple model, the averaged slow-subsystem is the switched 
system, 

-az 2 if zx < u sn 

*i = S d ' ( 13 ) 

-az 2 + j^pj if z\ > u sn 

-P(z 2 ~ Zi) if Zx < u sn 

^ = \ . ' ( 14 ) 



-P{z 2 - zx) + ^ if zx 



> Ms 



where, 



T(zx) = =(arctan( . c =) — arctanf . r )), (15) 

yfz\ +«sn V^l + Msn + «sn 

(see Figure [2] and Appendix |A]) . Now opposing contraction dynamics are explicit between 
terms on the RHS of Eq. [TH and Eq. [14] when - 



Zx > u s 



IV. CONDITIONS FOR MULTIRHYTHMICITY 

Non-monotonicity of C T is a necessary condition for multirhythmic bursting in circle/circle 
models because it means that P may support several isolated contraction mappings. Fluc- 
tuations in the contraction of P as measured by C T is the result of of near quantal spike 
addition in the active phase of bursting. Because the averaged system produces a monotonic 
C T , it cannot generate a saddle-node bifurcation of closed orbits (SCO) and therefore cannot 
support multiple coexisting stable solutions. 

Consider a closed orbit of the simple model 7 projected into the plane (ux, u 2 ) containing 
N spiking events. Dynamics on trajectories originating on the plane inside 7 is dominated 
by expansive spiking and they spiral outward. Dynamics on trajectories originating outside 
7 is dominated by contraction of the focus and they spiral inward. However, when initial 
condition (mx(0), m 2 (0)) lies some critical distance outside 7, the arc length of the resulting 
trajectory past u sn is long enough such that it contains iV + 1 spikes and its dynamics 
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may become net expansive, causing it to spiral outward, indicating a 'jump' in C T . If this 
occurs, this trajectory forms the inner boundary of the trapping region of a second oscillatory 
attractor. Since in Eqs. ( Fl3l H3J) quantal spiking is averaged over time, divisions between 
expansive and contractive annuluses about the equilibrium point are lost and P is a single 
contraction mapping for all initial conditions (figure EJ^c)). 

Concentric basins of attraction in the plane must be divided by an unstable invariant set. 
In Eqs. (0111 [5]), unstable periodic orbits are merely conceptual since spike addition occurs 
in a strictly discrete fashion. For a continuous system, unstable periodic orbits (UPO's) 
form separatixes for concentric basins of attraction. They are closed orbits that contain a 
dynamically unlikely attenuated action potential in the active phase. This can be seen in 
the biophysical models that support multirhythmic bursting^ 1 ^. UPO's in these systems 
correspond to the unstable fixed points of P which separate contraction mappings. 

C T is non-monotonic for Eqs. (El HI |5|) so moving toward some parameter sets causes 
its local minima and maxima to cross zero, indicating the creation or annihilation of a 
stable/unstable orbit pair via a SCO. The damping ratio of the slow subsystem is ( = 
{3/2y/a/3. As damping is reduced in the slow subsystem by (3 — > and a — > oo, dynamics 
move toward quasi-stable. This effectively reduces the average slope of C T pushing more of its 
extrema across zero causing SCO's to occur. Additionally, in circle/circle bursting systems, 
the firing rate past the SNIC scales as \Jb — b Bn where b is the bifurcation parameter and b sn 
is its value at the saddle node bifurcation (e.g Eq. [15] for the simple model). Since | 1 
increases as a — > 00 and — > 0, flow increases into u\ following these parameter limits. This 
allows spike addition to occur more readily as a function of initial condition (figure E]) • 

An interesting consequence of the mechanism underlying multirhythmic bursting in the 
simple model is that as a — > 00 and /3 — > and when d\ and o?2 are sizable, arbitrarily many 
stable bursting solutions can coexist. 

V. CLASSIFICATION OF POSSIBLE BURSTING BEHAVIORS 

From Section HVl it is apparent that a circle/circle bursting system can be reduced to a 
one-dimensional map that switches between two modes. This can be verified for complicated 
models of circle/circle bursting which produce highly nonlinear return maps, but preserve 
an alternating, sawtooth-like structure^. 
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FIG. 3. An idealization of P and C r for a circle/circle bursting model with a slow variable u. Stable 
and unstable fixed points are filled and open dots, respectively. As the slow subsystem becomes 
less dissipative and spiking events are capable of balancing contraction over larger areas of phase 
space. This effectively decreases the average slope of C T leading to the formation of multiple fixed 
points via SCO's (curved arrows). Parametric changes that drive trajectories of the slow subsystem 
further past u sn cause the number of spikes during an active phase to be increasingly sensitive to 
the initial condition. This effectively compresses C T increasing the number of separate contraction 
mappings via SCO's (straight arrows). P has three separate contraction mappings over the three 
sets of preimages within the dotten lines. 
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To explore the dynamics possible under this constraint, we introduce a simple piecewise 
linear map acting on a slow variable u that is depicted in figure [31 top. Starting with a local 
contraction mapping symmetric about the origin, the map alternates in both directions 
between two moduli, s\ and S2, which act over sets of size r\ and r^, respectively. Since 
the map defined by u n+ \ = s\u n is a contraction mapping, we confine | S2 j > 1 so that it 
is possible to have multiple basins of attraction. The average slope of the map is given by 
s ave = Sl l\+ s r 2 2 T2 ■ \s a ve\ < 1 requires the existence of at least one attractor under the map. 
Compared to P, the map is translated so that a central fixed point is located at the origin. 
Therefore, we consider the fixed point at the origin a nominal bursting solution with N 
spikes. Jumps to different contraction mappings indicate jumps to new bursting solutions 
with N zkk, k = 1, 2, 3, . . . spikes in the active phase. 

Now the limiting behavior of circle/circle bursting can be explored by inspecting the map 
for different values of s\ and S2- This analysis is carried between figure H] and table HI We 
note that the range of behaviors described by this simple system account for all the dynamics 
reported in previous multirhythmic bursting studies^ 1 ^ including the coexistence of chaotic 
attractors. However, those models are not topologically conjugate to our piecewise linear 
map since they can produce combinations of limit-sets seen in different parameter regimes 
of the one dimensional system under a single parameter set (e.g. the coexistence of a limit 
cycle and a strange attractor). 

VI. MULTIRHYTHMIC BURSTING IN BIOLOGICAL NEURONS 

We now provide evidence for the existence of multirhythmicity and multirhythmic burst- 
ing in invertebrate interneurons and neurosectretory cells. Cells R15, L3-L6, and L10 are 
identified neurons (neurons that are preserved animal to animal) located in the abdom- 
inal ganglion of Aplysia Californica. These neurons burst spontaneously in-situ with or 
without synaptic isolation. Additionally, they display the 'parabolic' bursting type as dis- 
tinguished by two key features. First, they are characterized by root scaling in firing rate 
during the active phase (figure E]) . Secondly, action potentials during bursting have after- 
hyperpolarizations that are lower than the voltage threshold for the active phase, ruling 
out a hysteresis mechanism for bursting 19 . Therefore we can conclude the circle/circle type 
mechanism for bursting well describes their dynamics. 
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TABLE I. The qualitative behavior of the piecewise linear map shown in figure H] is described for 
each parameter regime. The top section describes why some parameter values are irrelevant to our 
analysis. The bottom section describes relevant parameter ranges that can fully account for the 
dynamics witnessed in prior multirhythmic bursting models. 



Sl 


S2 


Dynamics 




\s 2 \ < 1 


S2 must have modulus greater than 1 


Si < 1 


s 2 < -1 


Non-biological because s ave < 


Sl > 1 


s 2 > 1 


Non-biological because s ave > 1 indicating that the system is a repeller 


< si < 1 


s 2 < -1 


Map is not well defined 


Si < -1 


S 2 > 1 


Chaotic and possibly multirhythmic since the full map is formed by a chain 
of sawtooth maps with only unstable fixed points 


-1< si < 


S 2 > 1 


Attractors are oscillatory fixed points and system is possibly multirhythmic 


< si < 1 


S 2 > 1 


Attractors are non-oscillatory fixed points and system is possibly 
multirhythmic 


Si > 1 


s 2 < -1 


Chaotic and possibly multirhythmic since the full map is formed by a chain 
of tent maps with only unstable fixed points 



R15 is the most heavily studied of the aforementioned neurons due to its particularly sta- 
ble activity and long burst period-. Numerous models^^^— have attempted to describe 
the mechanisms underlying bursting in this cell and many of these models are of the cir- 
cle/circle type^ 2 -. One class of models^ 2 ^ 2 ^ 2 ^ 1 ^ postulates that the slow variables underlying 
bursting are: are (1) the degree of voltage-gated activation of a slow, inward Ca 2+ current, 
Isi and (2) intracellular calcium concentration. Activation of l$i increases intracellular Ca 2+ 
and Isi displays voltage dependent activation and Ca 2+ dependent inactivation. 

We used a traditional 'current-clamp' technique to control applied current across the cell 
membrane while monitoring the membrane potential. Since l$i and intracellular Ca 2+ con- 
centration are directly and indirectly voltage dependent, current perturbations that change 
the membrane potential may influence these variables such that a multirhythmic cell is forced 
to switch attractors. Specifically, since our hypothesized mechanism of attractor switching 
involves spike addition, we aimed to change the number of action potentials in the active 
phase of bursting. 
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-l 1 

s, 

FIG. 4. Examples of the maps corresponding to the different parameter sets outlined in table HI 
Qualitative dynamics are classified by the following symbols: M - may support multirythmicity, 
C - chaotic, and NB — non-biological. Dashed sections of the maps show a single link in a chain 
of chaotic maps (sawtooth or tent map). The block for < s\ < 1, < — 1 is crossed out because 
the resulting map is not well defined. 

Models of R15 predict that outward current perturbations hyperpolarize the cell from 
its spiking threshold and when the cell is released from this hyperpolarized state it spikes 
vigorously ('rebounds'), leading to spike addition. Since activation of Isi is depolarization 
dependent, prolonged hyperpolarization completely deactivates Ca 2+ channels supporting 
Isi- Additionally, in the absence of an inward Ca 2+ flux, the intracellular Ca 2+ concentration 
decreases due to cytosolic buffering 2 - and the ion-channels supporting I$i are de-inactivated. 
When the inhibitory current is released the resulting depolarization quickly activates I$i 
allowing it to pass inward Ca 2+ current which maintains a high membrane potential causing 
spikes. Since intracellular Ca 2+ was initially very low, Isi must pass more Ca 2+ than usual 
before it is inactivated, resulting in a prolonged active phase. A depolarizing current has an 
opposing effect. Using 2 to 8 second, -0.5 to -2.5 nA hyperpolarizing current perturbations 
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5 10 15 20 25 30 35 40 

spike number 

FIG. 5. Early studies of R15 deemed the neuron a 'parabolic' burster since the spike-period profile 
resembles a parabola. However, it is now known that root scaling (grey line) better describes the 
frequency profile of this shape since after a SNIC the period of the resulting limit cycle scales as 
1/vA (equation [15] for the simple model). This figure confirms this behavior in an in-situ recording 
of R15. The instantaneous spike frequency during the active phase is shown with open circles. 

during the silent phase of bursting to add spikes to the subsequent active phase, we were able 
to show evidence for a multirhythmic behavior in four of eight bursting cells (see Appendix iBl 
for details of the experimental method). Figure [6] shows recordings in two different bursting 
neurons from separate animals that display a multirhythmic behavior along with the current 
perturbations used to induce an attractor switch. 

Section HVl predicts that an attractor switch in a multirhythmic burster will be character- 
ized by spike addition or deletion and an increase or decrease in burst period, respectively. 
[7J presents a comparison of slow- wave activity in [6(b) for two bursts prior and subsequent 
to the stimulus. Because coexisting bursting solutions are concentric in our models, the 
path-length of the solutions must differ resulting in distinct periods for each solution (see 
figure |2(a)). Note the existence of a constant phase advance developed after spike addition 
shown in figure [7(a) satisfies this prediction. In an attempt to elucidate the shape of the two 
coexisting attractors we plotted voltage versus its estimated derivative in figure UJb). Our 
multirhythmic model predicts that an unstable bursting solution containing an attenuated 
spike separates stable solutions. This separatix appears to exist, separating one attractor 
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from the other it proceeds for an additional spike in figure [7(c) . 

One may argue that the current perturbations we provide are too large both temporally an 
in terms of current amplitude to resemble a synaptic effect. However, the volume of charge 
moved across the membrane is a very superficial definition of synaptic efficacy. Smaller 
synaptic currents that rely on specific charge carriers (e.g. synaptic activation of calcium- 
specific channels) can be more effective in terms of altering specific neuron dynamics than 
our large chloride based currents. In the future, it would be interesting to repeat our 
experiment with direct perturbations to the hypothesized slow-subsystem via intracellular 
calcium uncaging. 

VII. POSSIBLE IMPLICATIONS OF MULTIRHYTHMIC BURSTING 

Since the time scales associated with short term memory, for example the multirhythmic 
motor memory suggested in£, are shorter than those associated with morphological synaptic 
plasticity, it is possible that neural systems employ of some activity dependent mult ist ability 
as a memory. Most proposed mechanisms of dynamic multistability in neural systems rely on 
some delayed feedback mechanism in a small neural circuit as the driving force in the creation 
of multiple coexisting attractors^— . Here we have shown that a sufficiently underdamped 
slow subsystem in a single circle/circle bursting neuron is enough to ensure the existence of 
multirhythmic behavior. 

We chose to show the existence of this phenomenon using invertebrate neurons because 
of their relative ease of experimental manipulation. However, as is put forth in sections 
IIVI and [Vj the dynamical mechanism underlying multirhythmic bursting is general and the 
multirhythmic regime occupies a non-finite range of parameter space; circle/circle bursting 
systems appear inherently suited for short term information storage. Circuits have been 
proposed that take advantage of this fact^. Even if a biological bursting neuron is not truly 
multistable, the shape of P indicates that the system will always have a non-monotoncity 
in the contraction of the vector field about the attractor (in C T ). This allows temporal 
amplification of perturbations to the system, even if it has only one true attractor (figure E]). 

In this report, we have proposed a basic dynamical mechanism for the existence of mul- 
tirhythmic bursting in the biophysical models that previously demonstrated this behavior. 
We then explored the range of dynamics possible in circle/circle type bursting systems. 
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Finally, we provided evidence for the existence of multirhythmic bursting in biological neu- 
rons. These experiments demonstrate a response to perturbations that is consistent with the 
dynamical model described herein. Further work concerning the nature of multirhythmic 
bursting neurons embedded within neural circuits is necessary to understand the importance 
of these findings in terms of short term memory. 
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Appendix A: Mathematical details of the simple circle/circle bursting model 
1. The quadratic integrate and fire neuron model 

Consider the quadratic integrate and fire model, Eq. [3j When u\ is treated as a parameter, 



where b Bn — I + U\. A saddle node bifurcation occurs when y^sn = indicating that the 

roots of the RHS of I A II have coalesced. When b sn < 0, these roots are real, v eq = ±\/h^, 
and are the equilibria of of IA1I . When b sn > 



which is an exact solution for the membrane potential v(t). 

From this we can calculate the periodicity of spiking once the neuron has moved to the 
tonic regime by finding the time for v(t) to go from v r to v c . Let t\ + to be the instant that 
a spike begins such that v(t\) = v r and t 2 + 1 be the instant that a spike terminates such 
that v(t 2 ) = v c . From|A2l 



v = b sn + v 2 , 



(Al) 





(A2) 




(A3) 
(A4) 




16 



Multistable Bursting 



Solving for t\ and t 2 , 




V I — 

ti = (arctan(^=)/ ^/b sn ) - t , 



(A5) 



(A6) 



Therefore, the inter-spike period (the period of a fast-oscillation) is given by T = t\ — t% 
which matches equation Eq. [15j 

2. Average bursting dynamics 

The slow subsystem of the simple model, Eqs. OH [5]), is T-periodically perturbed. Each 
time an action potential occurs in the fast subsystem, Eq. [3j there is an instantaneous 
modification of the slow variables. Let di represent the magnitude of this modification 
on a single slow variable, Ui, for each spike. Finding the average contribution of these 
modifications on the time derivative of Ui follows from the insertion of a delta function into 



The time constant, e, is dropped in this calculation because discrete changes to the slow 
variables are completely removed from contributions of time constants in the continuous 
dynamics defined byH]and[5j The validity of this result is confirmed by noting that, 



Therefore, when di/T{u\) is integrated over a single period, it produces an equivalent con- 
tinuous change in the slow variables as the discrete event does instantly per single period. 

Appendix B: Experimental materials and methods 

Aplysia californica were purchased from University of Miami Aplysia Resource Facility 
(Miami, FL) and maintained in a tank with filtered artificial sea water (ASW; Instant Ocean, 
Burlington, NC) at a temperature of about 19°C until used. Each animal was anesthetized 
prior to experimentation via injection about 50% of the animal's weight of isotonic MgC^: 
71.2 g MgCl 2 in 1 L of ASW solution containing (in mM) 460 NaCl, 10 KC1, 11 CaCl 2 , 30 
MgCl 2 , 25 MgS0 4 , and 10 HEPES (4-(2-hydroxyethyl)-l-piperazineethanesulfonic acid, pH 



Eq.ca 




(A7) 




(A8) 
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7.6). The animal was then pinned to a large dissection dish and opened rostrally to caudally 
along the dorsal midline. The abdominal ganglion was excised and pinned dorsal side up in 
a dish coated with Sylgard (Dow Corning, Midland, MI) in a saline solution of 30% isotonic 
MgCl2 and 70% ASW. The sheath of connective tissue covering the neurons of the ganglion 
was removed with fine scissors and fine forceps. The desheathed abdominal ganglion was 
then perfused with high-Mg 2+ , low-Ca 2+ saline, containing (in mM) 330 NaCl, 10 KC1, 
90 MgCl 2 , 20 MgSC-4, 2 CaCl 2 , and 10 HEPES (pH 7.6) to prevent synaptic transmission. 
Intracellular recordings were obtained using Clampex 8.2 software with an Axoclamp 2B 
amplifier (Axon Instruments, Foster City, CA, USA) in bridge mode using a microelectrode 
(10-20 MQ resistance) filled with 3 M potassium acetate. The membrane potential was 
amplified and digitized with a Digidata 1322A board (Axon Instruments) at a rate of 10 
kHz. Current pulses (-0.4 to -1.5 nA) of fixed length (2 to 8 seconds) were triggered externally 
using a square wave produced with an AFG3021 function generator (Tektronix, Beaverton, 
OR, USA) which itself was triggered by hand. 

We performed experiments on eleven bursting cells from eight animals. In order for a cell 
to be suitable for analysis, its activity needed to be stable and stationary so we could be 
confident that spike addition was experimentally induced and not the result of some intrinsic 
variability. Therefore, we required at least seven bursts to contain the same number of spikes 
before a stimulation was supplied and, if a mode switch occurred, that it was maintained 
for seven bursts subsequent to the perturbation. Eight of the eleven cells tested showed 
activity stationary enough to be analysed, and four of these eight cells showed instances of 
a sustained mode switch due to perturbation. In the four cells that showed no evidence of 
multirhythmicity, spike addition could only be maintained for the burst directly following 
the perturbation and the remaining bursts contained the nominal (pre-stimulus) number 
of spikes in the active phase. In cells that did show evidence of multirhythmicity, current 
pulses were not consistently able to induce a switch in bursting mode but our criteria for an 
attractor switch was met at least one time for each of the four cells, with three out the four 
cells showing multiple instances of a mode change. We saw no obvious qualitative correlation 
between characteristics of bursting (periodicity, duty cycle, number of spikes in the active 
phase, etc) and whether a cell was capable of a mode switch in response to perturbation due 
to the small sample size and simplicity of our study. 
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FIG. 6. Multirhythmicity in bursting neurons of Aplysia's abdominal ganglion, (a) A 'left-upper- 
quadrant' cell, most likely L3, showing bistability between bursting and tonic spiking. This be- 
havior is predicted in biophysical model o 15131 and was shown previously in R15^. Fundamentally, 
it is explained by the one-dimensional map P (section [V]) with N = 11 spikes per burst in the 
nominal (non-perturbed) mode and N = 1 spike per burst after the current pulse. Applied current 
pulses likely manipulate intracellular Ca 2+ concentrations, effectively moving the trajectory about 
in the plane of the slow variables, and allowing an attractor switch, (b) R15 showing bistability 
between two bursting solutions. Although the current pulse appears to have little effect, (c) shows 
the amount of spikes in the active phase of each burst. After the perturbation is applied, the 
number of spikes per burst increases by two, and then relaxes onto the new attractor containing 36 
spikes per burst. The active phase for pre- and post-perturbation is shown with the added spike 
highlighted. 
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FIG. 7. Attractor reconstruction from experimentally derived voltage traces, (a) Two bursts 
preceding and subsequent to stimulation were used for a comparative attractor reconstruction 
(bars in figure [6jb) denote the bursts used). A single burst period was defined by a the time 
between rising edges through -70 mV in the smoothed voltage trace. Voltage time series were low 
pass filtered using an exponentially weighted moving average with a width of 2 seconds to prevent 
spiking events from dominating the reconstructed attractor. Note the constant phase advance of 
one trace compared to the other, which is predicted by our model in figure [2j (b) Voltage plotted 
against its estimated derivative. This attractor shape can be directly compared to biophysical 
models of neurons which rely on a circle/circle bursting mechanism. A zoomed portion of the 
reconstruction is shown to highlight the existence of a separatix dividing the two attractors, an 
unstable closed orbit that contains an attenuated spike. The final portion of the active phase is 
thickened for each attractor to make this clear. 
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FIG. 8. 'Fading' multirhythmicity in L4 of Aplysia's abdominal ganglion, (a) Voltage trace and 
corresponding current perturbation. The three bursts following the perturbation contain an extra 
spike. This results from forcing the trajectory to pass through a bottleneck where dynamics 
are close to fixed which causes a lasting effect on bursting without the system actually being 
multirhythmic. (b) Fading multirhythmicity using the piecewise linear representation of P. Note 
that the map has areas that are close to forming fixed points. A trajectory can temporarily become 
caught here allowing a temporal amplification of a small perturbation (curved arrow), (c) Time 
series of the trajectory in (c). The ruins of an SCO allows the system to maintain N + 1 spikes 
per burst for a time before falling back into the true attractor with ./V spikes per burst. 
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